Dynamic polarization of graphene by moving external charges: random phase approximation 
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We evaluate the stopping and image forces on a charged particle moving parallel to a doped sheet of graphene 
by using the dielectric response formalism for graphene's 7r-electron bands in the random phase approximation 
(RPA). The forces are presented as functions of the particle speed and the particle distance for a broad range 
of charge-carrier densities in graphene. A detailed comparison with the results from a kinetic equation model 
reveal the importance of inter-band single-particle excitations in the RPA model for high particle speeds. We 
also consider the effects of a finite gap between graphene and a supporting substrate, as well as the effects of a 
finite damping rate that is included through the use of Mermin's procedure. The damping rate is estimated from 
a tentative comparison of the Mermin loss function with a HREELS experiment. In the limit of low particle 
speeds, several analytical results are obtained for the friction coefficient that show an intricate relationship 
between the charge-carrier density, the damping rate, and the particle distance, which may be relevant to surface 
processes and electrochemistry involving graphene. 
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I. INTRODUCTION 

The interactions of fast-moving charged particles with vari- 
ous carbon nanostructures have been studied extensively in the 
context of electron energy loss spectroscopy (EELS), typically 
using a scanning transmission electron microscope (STEM) 
with incident electron energies on the order of 100 keV. This 
technique has proven to be a powerful tool for investigating 
the dynamic response of carbon nanotubes 1,2 and, more re- 
cently, graphenei 3 - At such high incident electron energies, 
these studies have revealed important properties of both high- 
frequency a + tt plasmon excitations (~ 15 — 30 eV) and 
low-frequency tt plasmon excitations (~ 5 eV) in isolated 
single-wall carbon nanotubes 2 and in free-standing, undoped 
graphene. 3 While the studies on carbon nanotubes typically 
give plasmon dispersions at large wavenumbers (> 0.1 A -1 ) 
in the axial direction, 1,2 the study on graphene was performed 
with an electron momentum transfer close to zero, although it 
integrated over a significant in-plane component of the plas- 
mon wave vector. 3 In both cases, experimental data was found 
to be in good agreement with ab initio calculations. 2 ^ In ad- 
dition to ab initio calculations, methods employing an empir- 
ical dielectric tensor 5 and a two-fluid, two-dimensional (2D) 
hydrodynamic model for graphene 6,7 ° have also been able to 
reproduce the basic features of the a + tt and tt plasmon exci- 
tations in carbon nanostructures. 

For lower-energy external moving charges, recent progress 
has been made in measuring the dispersion of low-frequency 
plasmon excitations on solid surfaces using high-resolution 
reflection EELS (HREELS) with incident electron energies 
on the order of 10 eV. Such a measurement was performed 
on metallic surface-state electron bands, 8 and the results were 



interpreted theoretically using a dielectric -response model 
within the random phase approximation (RPA) that took into 
account the typically parabolic band structures of the sur- 
face states. 8 - 910 Furthermore, Liu et aZ.- ! have used HREELS 
to compare the low-frequency excitation spectra of doped 
graphene on a SiC substrate with the spectra of a metallic 
monolayer on a semiconducting Si substrate. At such low 
incident electron energies, the authors were able to measure 
the tt plasmon dispersion in a range of small wavenumbers 
(< 0.2 A -1 ) for a doped sheet of graphene with a high charge- 
carrier density. 1 1 This HREELS experiment, which provides 
the wavenumber-resolved spectra of low-frequency excita- 
tion modes in graphene with a high sensitivity to the doping 
levelr^ is more relevant to the parameter space in the present 
work than the STEM-EELS experiments. 

It is well appreciated that doping plays an immensely im- 
portant role in graphene's conducting properties, for which 
electron scattering on statically-screened charged impurities 
situated near graphene is one of the most important pro- 
cesses and is likely responsible for the famed minimum con- 
ductivity in undoped graphene. 12 ' 13 ' 14,15 ' 16 ' 17 ' 18 ' 19 In this con- 
text, important progress has been made in the development 
and use of the RPA dielectric function for low-energy exci- 
tations involving graphene's tt electron bands in the approxi- 
mation of linearized electron energy dispersion, which gives 
rise to the picture of massless Dirac fermions (MDF) i 20 ' 2 1,22,23 
This progress has opened up a range of interesting problems 
involving the interaction of graphene with external charges 
moving sufficiently slow that the MDF-RPA dielectric re- 
sponse theory can be applied. While an obvious application 
of this theory would be to interpret the HREELS experiment 1 1 
on graphene, another interesting application would be to the 
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study of slow, heavy particles moving near graphene. This 
latter application of graphene's MDF-RPA dielectric response 
theory would be relevant to studies of chemisorption of alkali- 
metal atoms, 24 friction of migrating atoms and molecule s 25,26 
moving near graphene, and ion transport in aqueous solutions 
adjacent to graphene when top-gating with an electrolyte is 
implemented. 27 Moreover, one could explore the application 
of low-energy, ion-surface scattering techniques 28 to graphene 
and other carbon nanostructures. There has also been re- 
cent interest in the directional effects of ion interactions with 
graphene-based materials, such as low-energy ion channel- 
ing through carbon nanotubes 29 and ion interactions with 
highly-oriented pyrolytic graphite, including implantation, 30 
channeling, 31 and ion-induced secondary electron emission 
from this target. 32 If applied to graphene, most of the scatter- 
ing configurations in these studies would involve impacts of 
slow, heavy particles under grazing angles of incidence, and 
many interesting parallels may be found with Winter's exper- 
iments on the grazing scattering of ions and atoms from solid 
surfaces^ 

We therefore wish to study the application of the MDF- 
RPA model to charged particles moving parallel to a single 
layer of supported graphene under gating conditions. In the 
wavenumber-frequency domain, (q, lo), the MDF-RPA model 
is applicable to graphene's polarization modes if the condi- 
tions q < 2k c and uj < 2e c /h are satisfied, where k c « a -1 is 
a high-momentum cut-off (with lattice constant a sa 2.46 A) 
and e c w 1 eV is a high-frequency cut-off validating the ap- 
proximation of linearized ir electron bands i 19 ' 21,23 For a point 
charge moving parallel to graphene at a fixed distance zq and 
constant speed v, the former condition will be satisfied only 
for distances zq > a, and hence we may neglect both the 
size of the particle and the size of the ir electron orbitals in 
graphene. The latter condition can be transformed into a re- 
striction on the particle speed by invoking the Bohr's adiabatic 
criterion and requiring that vJzq < 2e c /h. It is clear that with 
a gap on the order of 7 eV for graphene's a bands, particles 
moving at such slow speeds and large distances cannot excite 
the high-energy modes involving graphene's a electrons. 

Within the constraints of the MDF-RPA model, the main 
focus of this paper is on the stopping force and the dynamic 
image force acting on an external charged particle. We note 
that the stopping force is equal to the negative of the stop- 
ping power, which is defined as the energy loss of the exter- 
nal particle per unit length along its trajectory. 33 Meanwhile, 
the image force is a conservative force 34 that can strongly de- 
flect a particle's trajectory, especially for low particle speeds 
and/or small angles of incidence upon the target's surface. 33 
This was demonstrated not just for electron interactions with 
solid surfaces ^ 35 ' 36 ' 37 but also for ion 38 and molecule 39 grazing 
scattering from solid surfaces and ion 40,41 and molecule 42 ' 43 
channeling through carbon nanotubes. For example, in Ref. 38 
it was shown that both the stopping and image forces must 
be treated in a self-consistent manner in order to model ion 
trajectories and obtain ion energy losses that agree well with 
experiment results for the grazing scattering of slow, highly- 
charged ions on various surfaces^i 4 ^ A discussion of the stop- 
ping and image forces in the MDF-RPA model is therefore 



relevant to the current literature. 

In our previous work, 46 we have calculated the stopping and 
image forces on charged particles moving above graphene by 
assuming a high equilibrium density, n, of charge carriers in 
graphene and using a kinetic (Vlasov) equation to describe the 
response of graphene's tt bands within the linearized electron 
energy dispersion approximation!^ 4 ^ This semi-classical ki- 
netic equation (SCKE) model gave a relatively simple dielec- 
tric function for graphene that accurately described the ther- 
mal effect on plasmon dispersio n 47 ' 48 and allowed us to ana- 
lyze the contributions of plasmon excitations and low-energy 
intra-band single-particle excitations (SPEs) to the stopping 
and image forces. 46 However, it remained unclear how large 
the density must be to validate the semi-classical model and, 
more importantly, what effect the inter-band SPEs that lie 
beyond the capability of the SCKE model have. Therefore, 
the first goal of this paper is to determine the conditions un- 
der which the SCKE model is applicable at zero temperature 
by comparing the stopping and image forces obtained using 
the SCKE dielectric function 46 with those obtained using the 
MDF-RPA dielectric function . 21,22 ' 23 Furthermore, since we 
have found in Ref. 46 that a finite gap between graphene and 
the substrate strongly affects both forces in the SCKE model, 
the second goal of this paper is to examine the effect of a fi- 
nite gap in the MDF-RPA model. We note that the issue of a 
finite gap has become more important as increasingly diverse 
dielectric environments for graphene are studied^ 9 - 

Although we consider the MDF-RPA dielectric function to 
be a basic, parameter-free model that provides an adequate de- 
scription of the inter-band SPEs in graphene, the model never- 
theless has its shortcomings. For example, it ignores the local- 
field effects (LFE) due to electron-electron correlations^! and 
assigns an infinitely long lifetime to the electron excitations. 
The latter deficiency is often rectified in ab initio studies by 
applying a finite broadening, on the order of 0.5 eV, to the 
frequency domain for calculations of the loss function.— In a 
similar way, one can introduce a finite relaxation time, or de- 
cay (damping) rate, 7, to the MDF-RPA dielectric function 
for graphene using Mermin's procedure. 50,51 Since there are 
many scattering processes that can give rise to a finite lifetime 
of the excited tt electrons in graphene, an accurate determi- 
nation of 7 still presents a challenge. 51 ' 52 Therefore, the third 
goal of this paper is to treat 7 as an empirical parameter and 
investigate the effects of a finite damping rate on the stop- 
ping and image forces calculated with a MDF-RPA dielectric 
function modified by the Mermin procedure. This dielectric 
function, hereafter referred to as the Mermin dielectric func- 
tion, requires a careful extension of the MDF-RPA dielectric 
function derived for 7 = in Refs. 2122 to finite 7. Details of 
the Mermin dielectric function are given in Appendix A. 

The parameters of primary interest in this study are there- 
fore the equilibrium density of charge carriers in graphene, 
11, the graphene-substrate gap height, h, and the damping 
rate, 7. The equilibrium density is particularly important be- 
cause it determines the Fermi momentum of graphene's tt- 
electron band, kp — y/nn, and the corresponding Fermi en- 
ergy, ep = hvpkF, where vp « c/300 is the Fermi speed of 
the linearized tt band and c is the speed of light in free space. 
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In the case of intrinsic graphene (n = 0), the Fermi energy 
coincides with the Dirac point, £p = 0. In this paper, we 
consider a wide range of densities n > 0, expressed as a mul- 
tiple of the base value iiq = 10 11 cm~ 2 , under the conditions 
kp < k c and ef < £ c - 

The rest of the paper is outlined as follows. In section II, we 
present a theoretical derivation of the interaction of a general 
charge distribution with a layer of supported graphene. This 
derivation motivates the definition of the stopping and image 
forces for a point charge. In section EH, we compare the stop- 
ping and image forces in the MDF-RPA and SCKE models for 
the simple case of free graphene and a vanishing damping rate 
to determine the range of densities for which the SCKE model 
is valid. We then focus on the MDF-RPA model with a van- 
ishing damping rate, and in section IV we investigate the ef- 
fects of a finite gap between graphene and a SiC<2 substrate. In 
the simplified case of a zero gap, we provide analytic expres- 
sions for the stopping and image forces for intrinsic graphene 
and low particle speeds. Finally, in section V we consider the 
MDF-RPA model with a finite damping rate. After compar- 
ing the MDF-RPA model with experimental data to estimate 
the value of the damping rate, we consider the effects of the 
damping rate on the stopping and image forces with a special 
focus on the stopping force at low particle speeds. Note that 
we use Gaussian electrostatic units. 



II. BASIC THEORY 

We first give a brief generalization of the formalism de- 
veloped in Ref. 4 - to the case of a charge distribution with 
density p cx t(r, z, t). We assume that the charge distribution 
moves along a classical trajectory in a Cartesian system with 
graphene placed in the z = plane and with coordinates 
r = {x, y} in the graphene plane. In keeping with the re- 
flection geometry of ion-surface grazing scattering, 33 we as- 
sume that the external charge distribution remains localized 
in the region z > above graphene while a semi-infinite sub- 
strate occupies the region z < —h below graphene. Following 
Ref4&, we can express the induced potential $i„d(r, z, t) in 
the region above graphene using the Fourier transform (r — ► q 
and t — > u) as 

1 
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where e(q, oj) is the dielectric function of the combined 
graphene-substrate system and $[? xt (q, ui) is the Fourier trans- 
form of the external potential evaluated at the graphene plane, 
z = 0. The dielectric function of the system can be written as 

2ne 2 



e(q,w) = e bg (q) + 



■n(q,w), 



(2) 



where II(q, u>) is the polarization function for free graphene 
and £bg(<?) is the effective background dielectric function, 
which is expressed in terms of the substrate dielectric constant 



Esub as 
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We note that, instead of dielectric constant e su b, one may 
use a frequency dependent substrate bulk dielectric function, 
£sub(w), in order to include the effects of coupling between 
graphene's tt electrons and either the surface phonon modes 
in a strongly polar insulating substrate or the surface plasmon 
modes in a metallic substrate under the local approximation. 53 
In the former case, which includes a substrate with a single 
transverse optical (TO) phonon mode at frequency wto, one 
may use a dielectric function of the form^ 4 - 



e su b(w) = £oo + (e s - £oo) 
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u to - U ( UJ + *7to)' 
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where e s = e sub (0) and = lim^oo e sub (w) are the 
static and high-frequency dielectric constants of the substrate, 
respectively. In the latter case, which includes the high- 
frequency response of a metal, one may use the Drude dielec- 



tric function e su b(u;) 
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u 2 /[uj{uj + ijp)] with a pla 



frequency u> p and a damping rate j p . 

We limit the focus of this work to an insulating substrate 
in the static mode with a dielectric constant e SU b = e s , but we 
allow for an arbitrary gap h between graphene and the sub- 
strate. We note, however, that it is common in the literature to 
assume a zero gap j 16 ' 17,21 ' 22 for which Eq. (O gives an effec- 
tive background dielectric constant = (e s + l)/2. In this 
case, a simple description of the screening of electron-electron 
interactions in graphene can be quantified by the Wigner-Seitz 
radius, r s = e 2 / (e^Hvp)^- and free graphene can be recov- 
ered by setting e s = 1, and hence e° g = 1. Results provided 
for h — are therefore slightly more general than results pro- 
vided for h — > oo, which also characterizes free graphene by 
yielding eb g = 1 in Eq. (0. 

Next, we write <E>g Xt (q, u) = ^5(q, w), where the exter- 
nal charge structure factor 5(q, u>) is given by 
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S(q,«)= J dte lut J d 2 ve~^ r J dze~ qz p ext (r, 



z,t).(5) 



For a point charge Ze moving parallel to graphene with veloc- 
ity v and at a fixed distance z > 0, we find that 5(q, us) = 
2nZe 5(u> — q-v) e~ qz ° . In this case, the induced electric field 

E?" ld (r, z. t) = — V$?' ld (r, z, t) can be written as 



d z q (z-iq) e iq ' 



jq-(r-vi) c -q(z+z a ) 



E hHi( r ,M) = f^ 
- 1 



e(q,q- v) 
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where z is a unit vector perpendicular to graphene and q = 
<\/q. For stopping and image forces defined by F s = Ze v • 
E hid( r = v< > z = z o,t) andFi = Zez ■ E^ d (r=vt, z = z , t), 



respectively, where v = v/v, one obtain; 
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Note that we have used the symmetry properties of the MDF- 
RPA dielectric function e(q, ui) to simplify Eqs. (O and ©. 

For the comparison with the HREELS experiment in sec- 
tion V, we also define the total energy of the external charge 
reflected from graphene as 
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For a point charge Ze moving on a specular-reflection tra- 
jectory with p cx t(r, z, t) = Ze <5(r— V||t) S(z — v± \t\), where 
Vii and v± are the components of the particle velocity paral- 
lel and perpendicular to the graphene plane, respectively, the 
probability density for exciting the mode with frequency u> 
and wavevector q is^ 5 - 



{Ze) 2 v\q 



-1 



j-q-V||) +{qv±_' 
where we have set the reflection coefficient to unity 
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III. COMPARISON WITH SEMI-CLASSICAL MODEL 

In this section, we present the stopping and image forces 
calculated with dielectric functions from the SCKE model 46 
and the MDF-RPA model 21 ' 22 for free graphene (e° = 1) and 
a vanishing damping rate (7—5- 0). The results for both forces 
are normalized by Fq = Z 2 e 2 / (4zq), the magnitude of the 
classical image force on a static point charge a distance zq 
from a perfect conductor, to better reveal differences between 
the two models. 

Before proceeding, we note that the infinite upper limits 
of the q integrals in Eqs. dT) and © cause both forces to di- 
verge as the distance zq goes to zero. This behaviour, which 
also occurs in models of solid surfaces, 56 ' 57 should not be a 
major concern because the restriction z$ > a is necessary to 
ensure the validity of the MDF approximation. However, if 
one would like to extend the results for the stopping and im- 
age forces to include small distances, a standard procedure to 
eliminate the divergence at zq is to adopt a high-momentum 
cut-off. 56 ^ 5 - 7 . For studies of electronic processes in graphene 
with no external charges, it is common to impose a sharp cut- 
off at approximately fc e i 19,21 ' 23 However, there are other math- 
ematical methods for imposing a cut-off besides this sharp 
truncation of the q integration. 57 As discussed in Ref^, the 
use of an exponential cut-off function c~ 9/fec in Eqs. (0 and 
((H) would simply amount to a shift of the zq coordinate by 
a distance on the order of the lattice constant, which is some- 
times referred to as an effective image planed 5 - 7 - We therefore 



FIG. 1: The stopping force (a) and image force (b) normalized by 
F = Z 2 e 2 f{Azt) and shown as a function of the reduced speed 
v/vf of a proton (Z — 1) moving at a distance zq = 20 A above 
free graphene (h — » 00) for several values of the reduced charge- 
carrier density n/no, where no = 10 11 cm -2 . The thick and thin 
lines represent the results from the MDF-RPA and SCKE models 
with vanishing damping (7 = 0), respectively. 



evaluate the stopping and image forces using Eqs. (0 and ([8]) 
with infinite upper limits in the q integrals, and if one would 
like an estimate of the order of magnitude of these forces at 
small distances we note that a suitable shift of the zq axis may 
be chosen. It should also be mentioned that in an RPA model 
that takes into account the finite size of graphene's ir elec- 
tron orbitals, the divergence of these forces as zq — > can 
be removed by the resulting structure factor, which provides a 
physically motivated algebraic cut-off function.^! 5 - 7 - 

In Fig. 1, we compare the velocity dependence of the nor- 
malized stopping and image forces on a proton ( Z = 1) mov- 
ing at a distance zq = 20 A above free graphene in the MDF- 
RPA model (thick lines) and in the SCKE model (thin lines) 
for a broad range of densities. For intrinsic graphene (n = 0), 
note that both forces vanish in the SCKE model but they arise 
from inter-band SPEs in the MDF-RPA model. Furthermore, 
it can be seen that the results from the SCKE model agree with 
those from the MDF-RPA model only for high densities, and 
that this agreement is better for low particle speeds (v < vf) 
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FIG. 2: The stopping force (a) and image force (b) normalized by 
F = Z 2 e 2 /(4zo) and shown as a function of the distance zo of 
a proton (Z = 1) moving at a reduced speed v/vf = 0.5 above 
free graphene (h — ■> oo) for several values of the reduced charge- 
carrier density n/no, where no = 10 11 cm -2 . The thick and thin 
lines represent the results from the MDF-RPA and SCKE models 
with vanishing damping (7 = 0), respectively. 



than for high particle speeds (v > vp). The large difference 
between the MDF-RPA and SCKE models at high particle 
speeds is due to the presence of a plasmon line given by uo — 
uv(q) , 46 i 47 i 48 where u> p (q) = v F (q+q s )^ q/ {q + 2q s ) > qv F 
and q s = 4r s fcp is the Thomas-Fermi (TF) inverse screen- 
ing length. I3 ' 22 ' 46 Specifically, the energy loss for high par- 
ticle speeds in the SCKE model is dominated by the un- 
damped plasmon at frequency uo p {q), while the presence of 
the inter-band SPE continuum in the MDF-RPA model for 
ujjvp > max (q, 2kp — q) causes a strong Landau damping 
of the plasmon, 21,22 thereby producing much weaker velocity 
dependencies. However, even for these high particle speeds it 
appears that the SCKE model may be partially applicable un- 
der the condition z^kp ^> 1, which requires heavy doping of 
graphene in order to reduce the significance of the inter-band 
SPEs. 

The difference between the SCKE and MDF-RPA models 
at low particle speeds is further analyzed in Fig. 2. Using the 
same set of densities as in Fig. 1, we compare the reduced 



stopping and image forces on a proton (Z = 1) moving at 
a speed v — vf/2 above free graphene as a function of the 
particle distance, Zq. At such low speeds, one can see that 
the agreement between the SCKE and MDF-RPA models is 
better for the image force than it is for the stopping force. It 
follows from Fig. 2 that the condition z^kp > 1 may suffice 
as a rough criterion for the application of the SCKE model at 
low speeds (v < vp). This condition is far less restrictive than 
zokp ^> 1, which is required for the application of the SCKE 
model at high speeds (v > vp). We therefore discontinue 
further analysis of the SCKE model at high speeds, and turn 
our focus to analyzing various parameters in the MDF-RPA 
model alone. 



IV. MDF-RPA WITH VANISHING DAMPING 

In this section, we use the MDF-RPA dielectric function 
with a vanishing damping rate (7 — > 0) to evaluate the stop- 
ping and image forces. We first investigate the effects of a 
finite graphene-substrate gap, and then discuss two important 
cases with a zero gap: intrinsic graphene (n = 0) and vanish- 
ing particle speeds (v — > 0). 



A. Effects of a finite gap 

We now assume that graphene is supported by a SiC>2 sub- 
strate (e s w 3.9) and explore the effects of a variable gap 
height. It should be noted that the mean gap height between 
graphene and a SiC<2 substrate has been measured as 4.2 A^ 
which is comparable to the equilibrium distance of 3.6 A 
found in ab inito calculations between graphene and the top- 
most atomic plane of a SiC»2 substrate. 59 However, we note 
that h is defined in Ref. 46 as the position of an effective sub- 
strate surface plane where boundary conditions on the electro- 
static fields are imposed. Thus, while the measured and the- 
oretically obtained values of the graphene-substrate gap serve 
as a guide for the value of h, there is an uncertainty in h on 
the same order as the shift of the z axis discussed earlier. In 
this subsection, we consider the gap heights h — > 00 for free 
graphene, h = 4 A for a realistic value, and h = for the zero 
gap commonly considered in the literature. 

In Fig. 3, we compare the velocity dependence of the stop- 
ping and image forces on a proton moving at a distance zq = 
20 A above graphene for several gap heights and densities. 
For low particle speeds (v < vp), the gap height has a rela- 
tively small influence on the stopping and image forces that 
diminishes as the charge-carrier density increases and effec- 
tively screens out the graphene-substrate gap. The density, n, 
is therefore the most important parameter in the low-speed be- 
haviour of both forces. For higher particle speeds, the charge 
carriers in graphene are not as effective in screening out the 
graphene-substrate gap, and hence the gap height has a much 
stronger effect on the stopping and image forces. In particu- 
lar, Fig. 3 shows that for sufficiently high speeds (v vp) an 
increase in the gap height tends to increase the strength of the 
stopping force and decrease the strength of the image force, 
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FIG. 3: The stopping force (a) and image force (b) from the MDF- 
RPA model with vanishing damping (7 = 0) shown as a function of 
the reduced speed v/vf of a proton (Z — 1) moving at a distance zo 
= 20 A above graphene on a Si0 2 substrate (e s « 3.9). Results are 
shown for several values of the gap height h and several values of the 
reduced charge-carrier density n/no, where no = 10 11 cm -2 . 



FIG. 4: The stopping force (a) and image force (b) from the MDF- 
RPA model with vanishing damping (7 = 0) shown as a function 
of the distance zo of a proton (Z = 1) moving at a reduced speed 
v/vf = 6 above graphene on a Si02 substrate (e s w 3.9). Results 
are shown for several values of the gap height h and several values 
of the reduced charge-carrier density n/no, where no = 10 11 cm -2 . 



B. Intrinsic graphene with a zero gap 



but there is a range of moderate speeds for which this trend 
is reversed. Also note that for sufficiently high speeds, all 
MDF-RPA stopping and image forces approach the intrinsic 
case, n = 0. 



The effect of the gap height at high speeds is further ex- 
plored in Fig. 4, which shows the distance dependence of the 
stopping and image forces on a proton moving at a moder- 
ately high speed v = 6vp above graphene for the same gap 
heights and densities as in Fig. 3. In Fig. 4, it can be seen 
that the gap height has a strong influence on both forces for 
all distances. One may conclude that in the MDF-RPA model, 
as in the SCKE model, 46 any uncertainty or local variations in 
the gap height across graphene can lead to large fluctuations 
in the stopping and image forces, particularly for high particle 
speeds. 



We now take advantage of the simplicity of the MDF-RPA 
dielectric function for intrinsic graphene with a vanishing 
damping rat o 21 ' 22 to evaluate the stopping and image forces 
analytically for a zero gap. In this case, it is worth noting that 
the distance dependence of both forces can be factored out as 
Fq = Z 2 e 2 1 (4zg). For the stopping force, we find 



F 



Fo Ps 



1 



Pi 



-1/2 



(11) 



where p s = nr s /2 and v = v/vp. As seen from the 
thick solid curve in Fig. 1(a), this expression is subject to 
the velocity threshold constraint v > vp, which is a con- 
sequence of the inter-band SPEs yielding the loss function 
—^se(q,ui) > only for uj > qvpr^k For sufficiently high 
particle speeds (v ^S> Vp), one obtains an asymptotic form 
r i 2 

F s ° (tt/8) Ze 2 /(etz ) / (Hv) that is independent of 
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FIG. 5: The reduced image force J] from the MDF-RPA model 
with vanishing damping (7 = 0) for a proton (Z — 1) moving over 
intrinsic graphene (n — 0) on a substrate with dielectric constant 
e s and a zero gap (h — 0). The force is expressed as an effective 
background dielectric constant e£ g = 1/ [l + (42g/e 2 ) F°], and is 
shown as a function of the reduced proton speed v/vf and the actual 
background dielectric constant e[J a = (e s + l)/2. 



vp. It is interesting to note that the MDF-RPA stopping forces 
for all densities approach this high-speed asymptotic limit, as 
seen in Fig. 1(a) for free graphene and in Fig. 3(a) for various 
gap heights. 

The corresponding expression for the image force on a 
charged particle moving above intrinsic graphene, F®, is 
rather cumbersome. We therefore define the effective back- 
ground dielectric constant e£ g by writing the image force in the 

form F® = Fq (l/ebg ~~ l)' an d m Fig- 5 we present e£ g as a 
function of the particle speed v and the actual background di- 
electric constant e£ g — (e s + l)/2. We do this for background 
dielectric constants ranging from free graphene (e° = 1) to 
a HfC>2 substrate (e[J w 14). For vanishing particle speeds 
(v - 0), one finds e* B - e° g (1 + p.) « e a b 



% 22-, which 

2 Vp ' 



is a well-known result for the contribution of inter-band SPEs 
to the static limit of the MDF-RPA dielectric constant for in- 
trinsic graphene. 13 ' 22 For sufficiently fast particles (v ^> Vp), 
graphene becomes "transparent" and one recovers the case of 



a bare substrate, e 



b g 



e^g, with an accuracy to the order of 
e bg — e b° ~ v ~ 1 - Again, it is interesting to note that the MDF- 
RPA image forces for all densities eventually approach this 
high-speed asymptotic limit for intrinsic graphene, as seen in 
Figs. 1(b) and 3(b). 



C. Vanishing particle speed with a zero gap 

In this subsection, we consider the density dependence of 
the stopping and image forces from the MDF-RPA model in 
the limit of vanishing particle speed (v — > 0) for a zero gap. 
For sufficiently small speeds, the inset in Fig. 1(a) shows that 
the stopping force is proportional to the particle speed v, and 



hence we can define a friction coefficient rj through the equa- 
tion F s = —rjv. To evaluate rj, it is sufficient to expand the 
loss function 9f[— l/e(q, ui)] to the first order in uj and use the 
resulting expression in Eq. (0 to find the stopping forced We 
find the expression for the continuum of low-energy, intra- 
band SPEs, subject to the constraint q <2kp, to be 



2r, 



e bg e TF(9) 




(12) 



where £tf(<z) = 1 + q s /q is the TF dielectric 
function ! 13 i 22 ' 46 The friction coefficient is then given by 
rj = 2irhnZ 2 I(4:ZokF,r s ), where the function I(a,r s ) is 
defined in Eq. (4) of Ref.il 

The case z^kp <C 1 is particularly interesting because, un- 
like in the SCKE model, the constraint q < 2kp in Eq. ( fT2b 
causes I(4zokp, r s ) to remain finite even as z^kp — > 0. 46 
The friction coefficient for a charge moving very close to 
graphene in the MDF-RPA model is therefore given by 77 = 
2irhnZ 2 I(0,r s ), and is proportional to the charge-carrier 
density n. In the opposite case, z$kp ^> 1, one recovers the 
TF result for the friction coefficient given in Eq. (37) of Ref. 46 , 
which yields an asymptotic form rj = Z 2 h/(32zQy/nn) that 
is independent of vp. These two limiting cases for the fric- 
tion coefficient can be observed in Fig. 2(a), which shows the 
MDF-RPA stopping forces for a particle speed v = vp/2 that 
is very near the static limit. Recalling that the stopping force 
F s sa — \vpi~j is normalized by F = Z 2 c 2 /(4zq) in Fig. 2(a), 
note that all MDF-RPA curves are proportional to z 2 for short 
distances and fall off as z^ 1 for large distances. The transition 
between these two behaviours occurs around the peaks of the 
curves at z^kp ~ 1. 

In the limit of vanishing particle speed, the image force is 
closely related to the well-studied problem of static screen- 
ing of an external charge, for which the dielectric function in 
Eq. © reduces to e(q,0) = eL + — UJq), where TlJq) 
is the static MDF-RPA polarization function for graphene 
given in Appendix A i 13 ' 15 ' 21,22 Although the ui integration in 
Eq. is trivial with this dielectric function, the remain- 
ing q integration cannot be completed analytically. How- 
ever, expressions for the image force in the two limiting 
cases of z^kp can be given explicitly. For zo^f -C 1, the 
MDF-RPA image force reduces to the static limit for intrinsic 



graphene, F® = Fq 



/(£ 



1 



In the opposite 



ZL 2S. 

bg 1 2 v F j 

case, zokp ^> 1, one easily recovers the TF result for the 
static image force given in Eq. (38) of Ref. 46 , which yields 
asymptotically F\ ~ — Fq. Fig. 2(b) clearly shows the transi- 
tion between these two cases for the particle speed v = Vp/2, 
which is a good approximation to the static limit. Note that 
all normalized MDF-RPA image forces fall between the lim- 



its ( 1 + ff7>) -1 « -0.78 and - 
broad transition occurring at Zokp ~ 1. 



1 in Fig. 2(a), with a 
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V. MDF-RPA WITH FINITE DAMPING 

In this section, we use the MDF-RPA model with a finite 
damping rate 7 > (by means of the Mermin dielectric func- 
tion) to evaluate the stopping and image forces. Since the 
exact value of the damping rate is not known, we first treat it 
as a fitting parameter and come up with an estimate for 7 by 
comparing the MDF-RPA model with a finite damping rate to 
experimental data for the HREELS spectra of graphene on a 
SiC substrate. We then use this estimate to investigate the ef- 
fects of the damping rate on the stopping and image forces for 
free graphene (h — ► 00), with a special focus on its role in the 
friction of low-speed particles. 

It should be mentioned that introducing a finite 7 into the 
MDF-RPA dielectric function for graphene to create the Mer- 
min dielectric function, €m(q, 7)> is a non-trivial matter, as 
described in Appendix A. There is, however, a significant ad- 
vantage in using such a dielectric function in numerical calcu- 
lations of the stopping and image forces. Specifically, expres- 
sions for the MDF-RPA dielectric function in the strict 7 = 
case make numerical integrations of Eqs. ^} and © difficult 
as the boundary of the integration domain, < ui < qv, in- 
tersects the delta-like plasmon line. A small but finite 7, how- 
ever, broadens the plasmon line to allow for a simple numer- 
ical treatment of this delta-like behaviour. Throughout this 
paper, we have therefore reproduced the strict 7 = limit 
of the MDF-RPA dielectric function^i^ by using the value 
7 = 10~ 3 x vpko, where fco = ^/ttuo and uq = 10 11 cm~ 2 . 



A. Comparison with HREELS experiment 

To obtain a reasonable estimate for 7, we compare the 
MDF-RPA model with finite damping to the experimental 
data of Liu et al. 1 1 for the HREELS spectra of graphene on 
a SiC substrate. Since the focus of this paper is on slow- 
moving particles that are unable to excite graphene's a elec- 
trons, this HREELS experiment is more relevant to the present 
work than the EELS experiments on graphene and other car- 
bon nanostructuresi^i 3 . We note, however, that the effects of a 
SiC substrate on graphene are not as well understood as those 
of a SiC»2 substrate. Without entering the current debater 1 ^ 
we simply neglect any changes in graphene's 7r-band structure 
that may result from a hybridization of its it orbitals with the 
substrate and treat the gap height, h, as a free parameter. 

Since the measurements in ReL 11 indicate that the maxi- 
mum HREELS yields occur in the direction of near-specular 
reflection of slow incident electrons, it is appropriate to use 
Eq. ( fTOb for the probability density of exciting the mode 
(q, ui). However, since we are not aware at this time of spe- 
cific details of the experimental procedure that predominantly 
affect the low frequency range of the HREELS spectra via the 
prefactor in Eq. (fTOb, 11,55 we simply focus on the Mermin loss 
function 9f[— 1/ejvf (q, u>, 7)] and assume that it gives the cor- 
rect order of magnitude for spectral widths outside of this low 
frequency range. Therefore, in Fig. 6 we display a tentative 
comparison between the HREELS data 1 1 and the Mermin loss 
function with hj = 0, 200, and 400 meV and a gap height of 
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FIG. 6: The Mermin loss function (in arbitrary units) versus the en- 
ergy loss for graphene with a charge-carrier density n = 2 x 10 13 
cm -2 ) supported on a SiC substrate with static dielectric constant 
e s = 9.7 and a gap height h = 1 A. Model results are shown for 
damping rates hrf = 0, 200, and 400 meV, while symbols show the 
HREELS experimental data from Ref^i 



1 A for wavenumbers ranging from 0.008 to 0.102 A" 1 . The 
SiC subtrate is treated in the static mode with dielectric con- 
stant e s = 9.7, and the equilibrium density in graphene is set 
at n = 2 x 10 13 cm" 2 (hence e F « 570 meV and k F « 0.08 
A" 1 ) to match experimental conditions. Note that since the 
HREELS data is scaled arbitrarily, the Mermin loss functions 
for Try = 200 and 400 meV are scaled so that the maximum 
peak heights coincide with those from the experiment. For 
fvy = 0, however, the singular plasmon peak prevents such a 
scaling, and so the Mermin loss function is scaled by the same 
factor as for the Ivy — 400 meV loss function. 

In Fig. 6, the range q < kp is particularly interesting be- 
cause the Mermin loss function with h-f — exhibits three 
distinct features for these wavenumbers: a continuous spec- 
trum of intra-band SPEs for < ui < vpq, a continuous 
spectrum of inter-band SPEs for u> > vp(2kF — q), and a nar- 
row plasmon line at u> = 0J p (q) in the otherwise void interval 
vpq < u) < wf(2/cf ^ q)- The fact that these three features 
are not visible in the experimental HREELS spectra can be 
tentatively explained by assuming that a large enough damp- 
ing rate 7 exists, due to various scattering mechanisms, that a 
broadened plasmon line merges into the two regions of SPEs 
to form a single peak that follows approximately the original 
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FIG. 7: The peak positions of the Mermin loss function shown as a 
function of the wavenumber q/kF for graphene with a charge-carrier 
density n = 2 x 10 13 cm -2 supported on a SiC substrate with static 
dielectric constant e s = 9.7. Model results are shown for damping 
rates hrf = 0, 200, and 400 meV, as well as for gap heights h = and 
1 A. Filled circles show the HREELS experimental data from Ref." 



plasmon dispersion curve, oj = u> p (q). Note that a broadening 
of the plasmon line for the tvy = loss function does occur 
for the wavenumbers q = 0.077 and 0.102 A -1 as the plas- 
mon line crosses the boundary uo = vp(2kp — q) and enters 
into the region of inter-band SPEs, in which collective plasma 
oscillations decay into SPEs in a way that can be described 
by a finite Landau damping rate, 7^.— For the Mermin loss 
function with a phenomenological damping rate 7, however, 
a broadening of the plasmon line occurs for all q. Given that 
we have neglected the effects of temperature, the kinematic 
prefactor, and the reflection coefficient in Eq. ([Tol l, which all 
give rise to the low-frequency features in the HREELS spec- 
tra, Fig. 6 shows that a reasonably good qualitative agreement 
with the experiment can be achieved by using a gap height 
h = 1 A and a damping rate Jvy = 400 meV. 

In Fig. 7, we compare the peak positions of the model spec- 
tra in Fig. 6 with the experimental HREELS data. 11 In ad- 
dition to the best fit gap height h = 1 A, we also include 
the zero gap case to demonstrate the effect of the gap height. 
Note, however, that this effect is limited to shifting the spec- 
tral peak positions, and therefore does not significantly affect 
the estimate of the damping rate. In the long wavelength limit 
(q — » 0), it can be seen that the plasmon dispersion curve 
oj = u)p(q) for a vanishing damping rate exhibits the Jq be- 
haviour of a classical 2D electron gas, 21,22 but as the damping 
rate increases the plasmon dispersion curve falls off in a more 
quasi-acoustic manner. It is worth noting, however, that none 
of the curves are able to describe the trend of the experimen- 
tal data at long wavelengths. From Fig. 6, it can be seen that 
there is a strong coupling of graphene's it electron excitations 
with the non-dispersing Fuchs-Kliever surface phonon mode 
in the SiC substrate at the frequency w 1 16 meVj^^- We note 
that a dynamic treatment of the substrate phonon excitation 



through Eq. (O, the use of a g-dependent damping rate, and 
the inclusion of low-frequency features through the prefactor 
in Eq. ( fTob may be necessary for accurate modeling of the ex- 
perimental HREELS spectra in the full range of frequencies. 
However, further discussion on these points would go beyond 
the scope of the paper, and we merely conclude that a reason- 
able value for the damping rate hq is on the order of several 
hundred meV. 



B. Results for MDF-RPA with finite damping 

In this subsection, we discuss the effects of a finite damping 
rate on the stopping and image forces calculated with the Mer- 
min dielectric function for free graphene with several charge- 
carrier densities. Since it was found that a reasonable value 
for fvy is on the order of several hundred meV, we again con- 
sider the values h-f — 0, 200, and 400 meV, or equivalently 
7/ (vFko) ~ 0, 5, and 10, where fco = ^ttuq and no = 10 11 
cm~ 2 . 

In Fig. 8, we examine the velocity dependence of the stop- 
ping and image forces for the various damping rates and den- 
sities. For medium to high speeds (v > vp), it can be seen 
that the strength of both forces decreases significantly as the 
damping rate increases, and that this effect diminishes faster 
for the image force than for the stopping force as the particle 
speed increases. For low speeds (v < Vp), the fact that the im- 
age force is unaffected by the damping rate can be understood 
through Eq. ( IA5I) . from which it can be seen that the static 
limit of the Mermin polarization reduces to the static limit of 
the MDF-RPA polarization, U M (q, 0, 7) = U s (q). 

Fig. 8(a) also shows an interesting relationship between the 
damping rate, 7, and the density, n, at low speeds. One can see 
that the strength of the stopping force tends to increase with 
7 at a fixed density, but at a fixed 7 the stopping force tends 
to peak at some intermediate density. This trend is further 
confirmed in Fig. 9(a), which shows the stopping force as a 
function of the particle distance zq for low particle speeds. 
For high particle speeds, however, Fig. 9(b) shows that the 
distance dependence of the stopping force is affected more by 
the density than the damping rate. 



C. Friction coefficient 

The stopping forces in the inset of Fig. 8(a) suggest that, 
as for the case of vanishing damping, the concept of friction 
may be useful in the case of finite damping for low particle 
speeds. Again, we define the friction coefficient 77 through the 
equation F s = —rjV. Although we wish to focus on the case 
of free graphene (h — * 00), we provide analytic results for the 
slightly more general case of a zero gap and note that the case 
of free graphene can be recovered by taking e s = 1, and hence 
e bg = T° evaluate ?/, we expand the Mermin loss function 
for a zero gap to the first order in oj, which gives 
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FIG. 8: The stopping force (a) and image force (b) from the MDF- 
RPA model with a finite damping rate shown as a function of the 
reduced speed v/vf of a proton {Z — 1) moving at a distance zo = 
20 A above free graphene (h — > oo), Results are shown for several 
values of the reduced damping rate j/(vFko), where fco = y/nno, 
and for several values of the reduced charge-carrier density n/no, 
where no = 10 11 cm -2 . 



FIG. 9: The stopping force from the MDF-RPA model with a fi- 
nite damping rate shown as a function of the distance zo of a proton 
(Z — 1) moving at the low speed v — vf/2 (a) and the moderately 
high speed v = 6vf (b) above free graphene (h — > oo). Results 
are shown for several values of the reduced damping rate 7/(7>Ffco), 
where fco = y/nno, and for several values of the reduced charge- 
carrier density n/no, where no = 10 11 cm -2 . 



Eq. (13[ can then be substituted into the stopping force, Eq. 
0, to get an expression for the friction coefficient. Since the 
resulting integral cannot be evaluated analytically for an arbi- 
trary density n, we first consider the case of intrinsic graphene 
(n = 0), for which 

„2 



n(g, ii) 



1 



4fi 



(14) 



T 



(qv F ) 



and hence II S (q) = qj (4hvp). Using these expressions in Eq. 
( fT3] >, the friction coefficient for intrinsic graphene is given by 
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(15) 



where = 2zqj/vf, and Y v and H v are the Bessel function 
of the second kind and the Struve function, respectively. By 



combining the leading terms of the series expansions for large 
and small (o m this expression for rjo, one obtains the simple 
but surprisingly accurate formula (with a maximum error of 
approximately 3% at £ ~ 1) 
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(16) 



We now consider the friction coefficient r] for an arbitrary 
density n. In Fig. 10, we display the normalized friction co- 
efficient 77/770 for free graphene (eb g = 1) as a function of 
the reduced damping rate 7/(^^0)1 where fco = v/rmo, and 
the reduced charge-carrier density n/no, where no = 10 11 
cm~ 2 . It can be seen that the friction coefficient for a slow 
particle moving at a distance zq = 20 A above free graphene 
is a rather complicated function of 7/(7^^0) an d n/no, but 
the qualitative behaviour of this function can be understood. 
Recall from subsection III.C that the friction coefficient for 
a vanishing damping rate goes as 77 oc n for zofci? <C 1 and 
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FIG. 10: The friction coefficient r\ normalized by the friction coeffi- 
cient for intrinsic graphene 770 and shown as a function of the reduced 
damping rate 7/(«i?fco), where fco = ^/^n , and the reduced charge- 
carrier density n/no, where no = 10 11 cm -2 , for a proton (Z = 1) 
moving at a distance zq = 20A above free graphene (e° g = 1). 

as n oc 1/(zqu) for zokp 1, where kp — ^Jwn. For the 
distance zo = 20 A, the transition between these two lim- 
iting behaviours occurs at n/no ~ 80. For a finite damp- 
ing rate, we find in Eq. ( fToT ) that the friction coefficient goes 
as j] oc ^/(zqVp) for zoj/vp <C 1 and as rj oc Zq 2 for 
z ol/ v F 3> !• For me distance z = 20 A, the transition be- 
tween these two limiting behaviours occurs at ^/{vpko] ~ 9. 
Note how a saddle point develops at these two values for the 
density and damping rate in the surface plotted in Fig. 10. 



VI. CONCLUDING REMARKS 

We have presented an extensive analysis of the stopping 
and image forces on an external point charge moving parallel 
to a single layer of supported, doped graphene under condi- 
tions validating the massless Dirac fermion (MDF) represen- 
tation of graphene's 7r electron band excitations. These con- 
ditions require that the particle distance, zq, be greater than 
the lattice constant of graphene, and that the particle speed, 
V, satisfy Hv /zq < 2 eV. Calculations of the velocity and dis- 
tance dependencies of the stopping and image forces were per- 
formed within the random phase approximation (RPA) for a 
broad range of charge-carrier densities with three major goals: 
to compare the results with a semi-classical kinetic equation 
(SCKE) model, to examine of the effects of a finite graphene- 
substrate gap, and to explore the effects of a finite damping 
rate introduced using Mermin's procedure. 

With respect to the first goal, a comparison of the forces 
from the MDF-RPA and SCKE models in the regime of van- 
ishing damping has revealed that the latter model may be justi- 
fied for particle distances satisfying z^kp > 1 for v < vf and 
satisfying zqUf 3> 1 for v > vp. When combined through 
Bohr's adiabatic criterion, these conditions suggest that the 
SCKE model is valid only for heavily doped graphene with 



ep > Hv/zo, for which the effects of the inter-band single- 
particle excitations on the stopping and image forces are min- 
imized. With respect to the second goal, the effects of a fi- 
nite gap between graphene and a supporting substrate in the 
MDF-RPA model with vanishing damping have been found 
to be quite strong, particularly for medium to high particle 
speeds. These results have confirmed earlier findings from 
the SCKE model 46 and raise some concern over the com- 
mon practice of treating graphene with a zero gap when deal- 
ing with the dynamic -polarization forces on external moving 
charges. With respect to the third goal, we have made an effort 
to estimate the order of magnitude of the damping rate, 7, by 
providing a tentative fit of the MDF-RPA dielectric function 
modified by Mermin's procedure with recently published ex- 
perimental data for the HREELS spectra of graphene on a SiC 
substrate*^ To the best of our knowledge, this experiment is 
the only one that has been performed under the parameter con- 
straints of our model. With a suitable estimate of 7, we have 
found that both the stopping and image forces are primarily 
affected by the damping rate at moderate particle speeds. We 
have also shown that the combined effect of a finite gap and a 
finite damping rate can be important in modeling the plasmon 
peak positions of the Mermin loss function. 

In all calculations of the stopping and image forces within 
the MDF-RPA model, we have also found strong effects of 
the charge-carrier density, n, which are only weakened at suf- 
ficiently high particle speeds (v ^> vp). Although both forces 
were shown to approach the results for intrinsic graphene in 
this regime, it is likely that the diminished effect of graphene's 
doping level will be masked by excitations of graphene's a 
electrons as Hv/zq exceeds the a band gap (~ 7 eV). On 
the other hand, the parameter ranges considered in this pa- 
per are perfectly suitable for studying the friction of slow 
charges moving near graphene. The friction coefficient's de- 
pendence on the particle distance, the charge-carrier density, 
and the damping rate have been studied in detail for a zero 
gap, where analytic or semi-analytic results have been ob- 
tained. An intricate relationship has been found between these 
parameters and the friction coefficient, which may be of great 
practical interest for applications involving the concept of fric- 
tion, including surface processing and electrochemistry with 
graphene. 

There are several possible routes for extending the present 
work. Although not addressed in this paper, we have ex- 
amined the local field correction to the MDF-RPA dielectric 
function at the level of Hubbard approximation 63 and found no 
significant effects on the stopping and image forces. However, 
a recent treatment of the local field effects (LFE) using the 
GW method has led to significant improvements in the loss 
function of free, undoped graphene in the MDF-RPA model 
with vanishing damping. 4 Therefore, it would be desirable to 
treat the LFE using the GW method and explore its effects on 
the stopping and image forces. Furthermore, it would be de- 
sirable to extend the domain of applicability of the MDF-RPA 
model to small particle distances by including the finite size 
of graphene's ir orbitals. 20 Finally, the results of our tentative 
comparison of the Mermin loss function and the experimental 
HREELS spectra for graphene on a SiC substrate have opened 
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up the possibility of further modeling of this data in the low 
frequency regime. 
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APPENDIX A: THE MERMIN DIELECTRIC FUNCTION 

In this appendix, we extend the MDF-RPA dielectric func- 
tion for graphene given in Refs. l22l andl2ll to include a finite 
damping rate 7 > 0. Beginning with Eq. (3) of Hwang and 
Das Sarma, 22 we rewrite the polarization as H(q, lo + vy) — 
Il^(q,uj + i'j) + Il + (q,uj + ij), where it is assumed that q > 
and lo > 0. For a finite damping rate, the term IT (q, to + ij) 
is given by 



n (q,uj + %i) = 



g s gvqi 



where 



c 



lo + l"f 

VFq 



1. 



(Al) 



(A2) 



Note that the square root of a complex number z is defined 
as yfz = e2 lo s z , where logz = log|z| + iargz and the 
argument, arg z, is defined up to an integer multiple of 2tt. As 
with other multi-valued complex functions, the square root 
must be evaluated with respect to a single-valued branch, for 
which the value of arg z is selected from an interval of length 
2ir. For instance, the principal branch, denoted Arg z, selects 
the value of argz from the interval (— tt, tt]. In Eq. dAU , the 
square root must be evaluated with respect to the branch of 
arg z that selects values from (9 — 2ir, 9], where 9 = Arg £ is 
the principal argument of £. 
The term II 4 



n + (q,uj + ry) 



(q, lo + iy) is given by 
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(A3) 



where 



F(u) 



^Q{u 2 - 1) 



log (u VC + V((u 2 - 1)) (A4) 



and £ is defined in Eq. (IA2b . The square roots and logarithms 
in Eqs. (IA3t and (IA4b must also be evaluated with respect 
to the branch of argz that selects values from (9 ~ 2ir, 6], 
where 9 = Arg C Note that ^C(« 2 - 1) ^ V(Vu 2 - 1 
when dealing with specific branches of the square root, so Eq. 
( IA4l i cannot be simplified in the obvious manner. 



After introducing a finite damping rate 7, it is necessary 
to modify the polarization II(<7, lo + £y) using Mermin's pro- 
cedure to conserve the local number of charge carriers in 
graphenei^i The Mermin polarization function is then given 
by 



n M (g,w,7) = 



IL(q, cj + ij) 



1 - 



«7 



U(q,u + 17) 



n s (<?) 



(A5) 



where H s (q) = n(g, 0) is the static limit of the polarization 

Il(g, uj + ij) with u> — > and 7 — >• 0, given b y 21 i 22 



n s (<7) 



1 - 



LA'7. 



Q 
-ik F 



g s gvkF 
2tt Hvf 
if q < 2k F ' 



■SP- otherwise 

8k F 



(A6) 



Note that the piecewise definition of Eq. ( IA6b ensures that the 
square root and arcsine are real-valued. Using Eq. JA5t in Eq. 
(ffji, the Mermin dielectric function is then given by 

27re 2 

e M {q,u,y) = e bg (q) H 11^(^,0;, 7). (A7) 

q 

In the limit 7 — ► 0, it can be shown that Eq. JA7b is equiv- 
alent to the dielectric functions presented in Refs . 21 ' 22 t The 
difficulty in obtaining a compact form for the dielectric func- 
tion in the limit as 7 — > lies in reproducing the behaviour of 
the branch cut that naturally arises for 7 > — it is necessary 
to employ rather complicated piecewise-defined functions, as 
in Refs. 21,22 Alternatively, the compact Eq. iA7l may be used 
with a small, positive 7 to approximate this limit, or Eq. (IA71 > 
may be used with a realistic damping rate, 7, as originally in- 
tended. To this end, we describe how to implement the branch 
cut technique below. 

Most computer codes support basic arithmetic for complex 
numbers including exponentiation with base e, but only in- 
clude built-in functions for computing complex logarithms 
and square roots with respect to the principal branch. To com- 
pute logarithms and square roots with respect to the branch 
required in Eqs. (IA It . (IA3K and ( IA4b . it is necessary to imple- 
ment a function argbranch(z) that returns the argument of 
a complex number z in the range (9—2tt, 9], where 9 — Arg £. 
Using the built-in function at an 2 (y, x) , which returns the 
polar angle of the Cartesian point (x,y) in the range (— tt, tt], 
one may define 



Arg(z) : z 1 — ► atan2 (Im(z) , Rs ( z ) 



argbranch ( z ) 
Arg(z) - 2tt if Arg (z) >Arg(Q 
Arg ( z ) otherwise 



(A8) 



(A9) 



Since q > and lo > 0, Eq. ( |A9t correctly handles all values 
of z. The functions for the corresponding branches of the 
logarithm and the square root are then defined as 



logbranch ( z ) 



log ( I z I ) +iargbranch ( z ) , (A10) 
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sqrtbranch ( z ) : z i— » exp ( logbranch ( z ) /2 ) , (Al 1) 

where z | is the modulus of z, log ( ) and sqrt ( ) are 
the real-valued logarithm and square root, respectively, and 
exp ( ) is the complex exponentiation with base e. It is then 



a simple matter of evaluating Eqs. ( IAU . ( IA3I ). and ( lA4b with 
respect to these branches of the logarithm and square root. 
Note that argbranch (() must return Arg (0 to yield the 
correct result. 
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